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ABSXEWCC,,. 


A one^dimensian^X,. electron' .fluid-, model., is- investigated using 
the mathematical- methods of modern fluid turbulence theory.,,, lion- 
dissipative equilibrium, . canonical distributions, are determined in a- 


phase , space, whose ^co-ordinates , e ,jth,e ^eal, ..aftd .q’^agiripry.y.partvr., pf , 


the Courier,.. coefficients for the field* variables. • Spectral densities 


are calculated,,, .yielding, a >. wavenumber., .electric,, field, energy spectrum 

-2 

proportional .to* to ’for largo wavenumbers. The equations, of motion,, 
are numerically integrated .and the resulting spectra, are found .to 
compare well with the theoretical predictions. 



I. INTRODUCTION 


Though not yet a completed program, the statistical theory of 
fluid turbulence has made remarkable strides in the last two decades. 
Many of these advances hare not been digested by the plasma turbulence 
community, however. Plasma turbulence theory, consequently, has 
developed in a considerably less convincing conceptual framework and 
has achieved less detailed verification of its predictions in experi- 
ment and numerical simulation than fluid turbulence has. 

It Is of interest to see how many turbulent plasma situations 
can be treated by methods borrowed from fluid turbulence theory. 

Though no survey of the plasma turbulence literature is attempted 
here, we remark upon some recent treatments of incompressible magneto- 
hydrodynamic turbulence 1-7 in both two and three dimensions. The 
present article concerns a compressible situation: the electrostatic 

turbulence of an electron fluid plasma model. The challenging problem 
of the Vlasov plasma, despite its enormous literature, does not appear 
to us to be quite ripe for an approach through modern turbulence 
techniques, except possibly in the so-called "weak turbulence" limit, 
which has been extensively studied. 

The principal concern in fluid turbulence has been with high 
Reynolds number situations. . For these, the nonlinear terms are so 
much larger than the linear ones that the dominant physical effect 
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can be accurately said to be the transfer of the excitations from 
one spatial scale to another. Often, migrations of energy across 
orders of magnitude in spa! i al scale are involved. These migrations 
can "be toward either shorter or longer wavelengths, depending upon 
the circumstances. There is an accompanying emphasis in the theory 
on the wavenumber spectrum, on energy transfer functions, on cascades 
through wavenumber space, and so on. An unanswered question (almost 
an unasked one) in electrostatic plasma turbulence concerns the 
extent to which the phenomena also involve large migrations of energy 
from one spatial scale to another. With few exceptions, theories 
have been formulated in ways that implicitly assume that no such 
migrations occur, and provide no mathematical machinery to describe 
the transfer, numerical simulations typically have not involved 
sufficiently fine spatial resolution to see such migrations, and the 
intensity of the search for a stable, quiescent equilibrium has often 
led simulators to stop computing as soon as spatial fine structure 
has begun to develop. Spectral-method computing techniques , 8 ” 11 
which are well suited to computing simultaneously on a wide range of 
spatial scales in Navier-Stokes and magnetohydrodynamic fluids, have 
not yet had much impact on plasma simulation. Configuration-space 
computations are less than satisfactory for following evolving 
turbulent fields. 

There are some good reasons for the neglect. Plasmas do not 
obey the same dynamical equations over as many orders of magnitude 
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in spatial scale as Kavier-Stokes fluids or 'stellar- interior MID 
fluids typically do. For example, while MHD fluid descriptions may 
satisfactorily represent the largest-scale dynamics of a tokamak 
plasma, on the scale of the ion gyroradius there are important 
phenomena which are not contained in MHD at all. The incorporation 
of quite different mathematical descriptions within the calculation 
of a single turbulent situation lies out of reach of known procedures. 
All this implies that most situations about which we can make logically 
compelling quantitative assertions will he "unrealistic” in the sense 
that they will omit phenomena of importance for real experiments. 

One of two alternatives must be chosen: (l) continue to ignore the 

possibility that interesting plasmas are, like other fluids, turbulent 
in a sense that involves transfer of the energy in the field from 
one, spatial scale to another; or (2) study highly simplified models 
that, while certain never to tell the whole story about a turbulent 
plasma, may at least tell a part of it accurately. In the present 
article we opt for the second choice as being historically the more 
likely one along which progress may be expected. 

In Sec. II we present the details of the model of an electron 
plasma to be studied. It is important to recognize that in restricting 
attention to an electron plasma with uniform immobile ions, we will 
be omitting one of the more popular phenomena of the day: "Langmuir 

solitons", and other related situations in which the ponderomotive 
force is involved in a fundamental way. In Sec. Ill the statistical- 
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mechanical predictions are derived. Section IV is a presentation of 
some preliminary numerical tests of the material presented .'.n Sec. 
III. Section V summarizes and discusses the results. 
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II. THE MODEL 


We take the simplest possible model in which 'electrostatic 
effects in a plasma will be exhibited: the one-dimensional electron 

plasma with a locally adiabatic equation of state. If the electron ■ 
density is n, the fluid velocity is v, the electric field is E, and 
the ratio of specific heats y is 2, we have, in dimensionless units, 

H + & (nv) = 0 » (1) 
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( 2 ) 


SE 

Sx 


= 1 - n 


Vie may replace Poisson's equation, Eq. (3), by 


8E 

St 


nv 


(5) 


<4) 


•which is an equation for the time advancement for the electric field, 

/ 

if we invoke Poisson's equation as an initial condition. Equations 
(1), ( 2 ), and (4) will then preserve the Poisson relation in time. 
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n, E, and v in Eqs. (l) — ( 4 ) are functions of only one space 
coordinate x, and the time t. n is measured in units of n , the 
equilibrium average number density of the electrons and (immobile) 
positive ion background. Times are measured iri units of the inverse 
electron plasma frequency o' and velocities in units of the electron 
thermal velocity; lengths are in units of the Debye length. The 
mechanical pressure gradient enters in the second term on the right- 
hand side of Eq. (2). Any relation between the pressure and density 
other than (p/n ) = constant would have complicated the pressure 
gradient term in Eq. (2), but we believe that no qualitative restric- 
tions are introduced in the physics by this particular choice of y - 2. 
The term v is a dimensionless viscosity which can be thought of as 
the reciprocal of a large Reynolds number. The limit v = 0 will be 
called the "non-dissipative" limit, and is the only limit in which 
Eqs. (l) — ( 4 ) are a conservative system. [A linear term proportional 
to v and containing a collision frequency is also possible in Eq. (2).] 
There are typically two standard situations in which the pos- 
sibility of a tractable steady-state statistical theory of the Eavier- 
Stokes equation has existed; (l) the non-dissipative initial value 
problem; and (2) the driven, dissipative, steady-state problem in 
which some external agency supplies energy at a given spatial scale 
while dissipation occurs at some shorter spatial scale as a consequence 
of viscosity. The second problem is much more physically realistic 
than the first, but valuable insights are obtained from studying the 
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non-dissipative problem as a preliminai'y. If v:e want to study the 
driven, dissipative case for Eqs. (l) — (4), it is necessary to add 
to the right-hand side of Eq. (2) a driving term which represents an 
external source of the excitations, and to allow v to be non-zero. 

The really new feature present in Eqs. (l) — (4) which to our 
knowledge has not yet been incorporated into either MHD or Ifevier- 
Stokes theories is the finite compressibility. Incompressibility has 
been a useful simplifying feature of PTavier-Stokes and MID fluids as 
far as turbulence goes, but its introduction into a one-dimensional 
plasma model would render the model meaningless. Compression is 
essential for electrostatic fields to fluctuate, since they are 
generated by density inhomogeneities. 

That some qualitatively new and unanticipated features are to 
be expected as a consequence of the finite compressibility can be 
seen by looking at the right-hand side of Eq. (2) when v - 0. Compare 
the two surviving terms with respect to order of magnitude, for an 
excitation of typical wavenumber k. For a given perturbation of the 
density of magnitude say, Poisson's equation tells us that the 

magnitude of the electric field to be expected is ~ fih^/k. The 
magnitude of the pressure gradient term, on the other hand, will be 
~ kSn^.. Thus at very long wavelengths (» a Debye length) the pres- 
sure gradient term will be negligible, while at very short wavelengths 
(« a Debye length) the electric field will be negligible. The 
excitations will go from being very like a nonlinear cold plasma 
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oscillation at long wavelengths to very like a nonlinear sound wave 
at short wavelengths, with a continuous transition in the spatial 
scales of the order of a Debye length. (This is one respect where 
the behavior of the model will differ sharply from that of the Vlasov 
plasma, which, lacking collisions, will also lack the coherence at 
short "wavelengths to exhibit the sound waves. ) 

The short wavelength behavior of Eq. (2), then, should be that 
of the Euler equations of a compressible perfect gas. Compressible 
turbulence usually occurs in gases as a perturbation on incompres- 
sible turbulence, and has been studied in isolation very little. One 

/ 

thing our experience with the Euler equations of compressible flow 
strongly suggests as a possibility, however, is a tendency toward 
the formation of shocks, or discontinuities. A picture is suggested 
in which long-wavelength turbulent plasma oscillations feed short- 
wavelength sound waves, which steepen towards the development of 
shocks, which in turn are eaten away by the action of the dissipation. 
The shocks may be very weak ones, and it is imaginable that dissipa- 
tion might destroy them before their formation becomes complete. 

This is a gross, highly conjectural picture for the short wave- 
length behavior, but is reminiscent of that which is known to occur 
for Burgers’ equation. 18 " 14 [Burgers' equation is just Eq. (2) minus 
the first two terms on the right-hand side. ] Our program here is to 
first apply the statistical methods of turbulence theory to the 
non- dissipative limit of Eqs. (l) — ( 4 ), and then to test some of our 
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predictions by numerical solutions. ' Inclusion of finite dissipation 
is a larger matter, and is deferred to a later publication. 



III. 


STATISTICAL FOFiMULATIOH 


We impose periodic boundary conditions over a length L and 
Fourier-decompose, in a plane-wave expansion, v, n, and E: 

v - \ v(k* t) exp (ikx) , 

n = n(k, t) exp (ikx) , 

E = E(k, t) exp (ilex) (5) 

« 

where is over all wave numbers 2nm/L ; with m = 0, ±1, ±2, ... It 
is e'ssential for complete generality to include m = 0. It can be 
shown that n(k - 0, t ) = constant = 1, but it is not true that we 
can automatically set E(k = 0, t) = 0 or v(k - 0, t) = 0. It is a 
common belief that the option of setting the k = 0 components of E 
and v identically zero is the only option available or the most 
natural and acceptable one. Appendix A is devoted to clarifying and 
correcting this widespread misconception. 

The equations for the Fourier coefficients are (with time 
arguments suppressed for economy) 

— ■ ip- + ik £ n (r) v(p) = 0 > (6) . 

p+r=k 
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pY ffi - + E_ v (p) = - e oo - ikn oo 


( 7 ) 


p+iv-k 


ikE(k) = -n(k) 


(8) 


Equations ( 6 ) and (7) hold for all k. Equation ( 8 ), however, holds 
only for k ^ 0. The use of Eq. ( 8 ) will eliminate' E(k) in favor of 
n(k) for k / 0, hut an equation is needed E(0). From Eq. (4) 


afS* = E n(k) v(-k) 


( 9 ) 




Equations ( 6 ) — ( 9 ) are, in general, an infinite set of ordinary dif- 
ferential equations for the Fourier coefficients in the infinite 
Fourier series and as such are intractable. If further progress is 
to be made, some simplification has to he made. The simplification 
made in a linearized theory would he to replace the convolution sums 
in Eqs. ( 6 ) and (7) with the two terms where either p = 0 or r = 0, 
thus eliminating all Fourier coefficients from the equations except 
those with wavenumbers k and 0. The resulting equations are then 
easily solved. Basically the assumption in the linearized theory is 
that all non-zero wavenumber Fourier coefficients are small compared 
to the zero wavenumber Fourier coefficient n(0). The very nature of 
strong turbulence, however, does not allow one to make this approxi- 
mation. The standard practice in fluid turbulence is rather to 
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truncate the Fourier series at a large, but finite, "wavenumber Y • 
This truncation can be reasonably justified for the dissipative 
problem. In that situation one envisions taking k so large that 
for any -wavenumbers larger than k , dissipation wipes out any 

niaX 

Fourier coefficient before it attains any significant value. This 
is not true, however, in the non-dissipative limit. Nevertheless 
insights into how the .nonlinear ..terms ..act to ..shuffle energy between 
the various length scales can be obtained through this truncation 
just as these insights are gained in 'the Navier-Stokes and MHD 
problems . 

The difficulties encountered in trying to solve the large, but 
finite, truncated system of equations is not significantly different 
from those encountered in the dynamics for a large, but finite, 
system of point particles (also described by ordinary differential 
equations) in statistical mechanics. In a sense, what we intend to 
do here is a statistical mechanics of Fourier coefficients. 15 In a 
phase space consisting of the real and imaginary parts of independent 
Fourier coefficients we can write down a Liouville equation and ask 
for a stationary solution. As we know from elementary statistical 
mechanics, such solutions depend critically upon the "constants of 
the motion." We have discovered three conserved quantities for Eqs. 
(6) — ( 9 ). Two of them are the energy £. and the current (or momentum) 
J, which in terms of the Fourier coefficients are 
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& = h £ 6(k, + k 0 + k,) n(> ) v(k 0 ) v(k,^ 

• * k x k 2 k^O 1-3123 

+ |E |v(lOj 2 +§E |n(k)i 2 +|C |E(k)l 2 , (10) 

k^O k/0 k/0 

J = Z! n(k) v(-k) . (11) 

k/0 

Hotice that the total current (momentum) is not conserved in time, 
hut only the current associated with the k ^ 0 modes. 6 (k^ + k^ + k^) 
in the first term of & is a Kronecker gelta which has the value 1 
whenever k^ + kg + k^ = 0-and the value 0 whenever + k^ + k^ ^ 0. 

As a result the first sum in £ is a sum over all sets of three wave- 
numbers with 0 < Ik. I S k , i = 1, 2, 3 allowed by the periodic 
boundary conditions which add to zero. This first term we mil 
colloquially call the "triples". The remaining sums in £ and the 
only sum in J are over all wavenumbers allowed by the periodic 
boundary conditions with 0 < [k| s k^ y . To our knowledge £ is the 
first non-quadratic constant of the motion for a truncated Fourier 
system to appear in a turbulence calculation. (See Appendix B. ) 

The constancy of the current reduces the equations for the 
k - 0 Fourier coefficients to a pair of linear equations 


- - E <°> - 


( 12 ) 
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^i=.v(0)+j (13) 

where n(0) = 1 has been used. As a result 

Q. = CE(0)] 2 + [v'(0) +#] v(0) (14) 

is also a conserved quantity. Equations (12) and (.13) have the exact 
solutions 

v(0) = -J + (v Q + J) cos t - E q sin t , (15) 

E(0) = (v + J) sin t + E q cos t (l6) 

where v Q and E q are the initial values of v(0) and E(0), respectively. 

We believe that £, J, and Q, are the only conserved quantities 
which survive the truncation of Eqs. (5)~(9) to a large, hut finite, 

t 

number of Fourier coefficients. A condition that we cannot show and 
do not believe is true is that the number density n(x, t) formed from 
the finite Fourier series will remain non-negative at all configura- 
tion space points for all time even though the initial number density 
was non-negative at all points. It would not be totally surprising 
that a truncated Fourier series would yield a non-physical number 
density after a finite time in the non-dissipative limit. 
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The quantities -which receive the roost attention in turbulence 
calculations are' the wavenumber energy spectra. These spectra here 
can be calculated as moments of a canonical ensemble constructed in 
the usual manner from the conserved quantities with an inverse temp- 
erature (Lagrange multiplier) for each conserved quantity. Unfor- 
tunately the triples male the calculation of moments using the full- 
canonical distribution impossible. Vie, therefore, consider a "weak 
turbulence" limit in which the product of three Fourier coefficients 
can be considered small compared to products of tiro coefficients. If 
& -represents g without the triples and Cl = [e(0)j + [v(o)] then 

we propose to use the canonical distribution 

I = T| exp (-G£ - pj - rQ*) (17) 

to calculate the moments of Fourier coefficients. T] is a normalizing 
constant determined by the normalization 

-J’ D dX = 1 (18) 

where j* dX is an integral over all the independent real and imaginary 

# -X- 

parts of Fourier coefficients. Since n(k) = n (-k) and v(k) = v (-k), 
independent Fourier coefficients are associated with only half the 
wavenumbers. Furthermore E(k) and n(k) for k ^ 0 are not independent 
since they are related through Poisson's equation, a, f3, and y play 
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the role of inverse temperatures and are constained by the requirement 
that D be normalizable. 

•)*- 

We can make a further simplification. Q involves only the 

* ■ , 
k - 0 inodes while E and J involve only modes with k f 0. Since we 

already have exact solutions for the k = 0 modes and are interested 

in expectation values fo'r the energy in the k ^ 0 modes we can use, 

for the purposes of the calculation 

D = 71 exp (-QS* - pj) ( 19 ) 

* 

Furthermore, since £ and J are sums of terms indexed by a single 
wavenumber, D factors into a product of distributions , one for each 
k. The single -1: distribution is 


f k tn r (k), n i (k), v,(k), v^k)] 


= \ exp 


-a|v^(k) + v?(k) + Jl + 


[u^(k) + n?(k)j 


^ [n r (k) v r (k) + n i (k) v i (k)] 


( 20 ) 


where n(k) = n r (k) + in^k) and v(k) = v r (k) + iv t (k). The normalizing 


constants are given by 
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a 

2 

TT 


a 1 + 




a 


( 21 ) 


The requirements that the 


distributions be normalizable are that 


Ct> 0 , 


a 


1 




for all k. 


( 22 ) 


The latter condition can be met for all k if it can be satisfied for 
k 

max 

Expectation values can now be computed. 


<*V(k)> - <nf(k)> 


1 


i 


2 

a 

+ jli 

1 eil 

| a 


» ' 

l k j 



( 25 ) 


<v^(k)> = <v 2 (k)> 



( 24 ) 


< n r ( k ) \( k )> = <n ± ( k ) v i (k)> = 


P 1 


2a 

tafl + -^1 

P 2 


L 1 k 2 i 

a 


( 25 ) 
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Hence the wavenumber energy spectra are Given by 


<l' r (k) f) 


i + p; 

ry 2 
a 


a 


1 + 


k‘ 




,,2 


u 


<|n(k) f) - 


a | 

[l + -~j 

_ ef 1 


i k 2 / 

a 


and 


<|E(k)j 2 > = 


X 


-> 

aj 

[l + -ij 1 

p 2 " 

j <2 j 

i " ~cc 


l k J 



( 2 ^) 


( 27 ) 


( 23 ) 


a and (5 are determined by the condition that the expectation values 
* 

of £ and J match their initial values; that is, a and & are solutions 
of the algebraic equations 


* 

e 



i 

a 


, i °q + fc a ) + (^r)fc g 

v ; lai J 


J = 


. £ § k 2 

^ 0 “an + K 2 ) - (£)k 2 


( 29 ) 


( 30 ) 
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The electric field energy spectrum, as a function of vavereS: 
has one basic shape: concave down, proportional to k~‘ for kX- I 

and approaching the constant valnc a ^ for It « 1 vith a contir. 
transition region in between. 

In the particular case that J = 0 we have that p = 0 and. the 
spectra reduce to the simpler forms 

<ls(k)| 2 > OD 

a(l + k ) 

and 

<|v(k)| 2 >»i . (32) 

It is this particular case which we will examine numerically in the 
following section. 



Equations ( 6 ) ~ ( 9 ) are solved numerically ns its 3 the spectral 

method of Orszag and rattorson .''""' 1 G The essence of the spectral 

method is to evaluate the convolution sums "which appear in these 

'-equations by using the convolution theorem and a fast Fourier 

Transform. Unfortunately the convolution theorem does not recognise 

truncations in k-space. Modes v, T ith Ikj > k can be generated, 

max 

leading to what are termed aliasing errors. Fe have eliminated 

these aliasing errors hy setting to sero the outer one- third of our 

computational e.rray at the end of each time step. That is, if IT is 

the total number of configuration space grid points, we have retained 

only the 2 K /3 wavenumbers of smallest magnitude out of the possible 

N Fourier modes by setting k = (lT/3 ) k . where k . is the smallest 

non- zero wavenumber allowed by the periodic boundary conditions. 

The time advancement is done using a fourth-order Adams - 

Bashf or th -Moult on. predictor- corrector with the first two time steps 

being calculated by a fourth-order Bunge -Kutta method . 12 

Table I gives a list of parameters for three typical runs in 

which there is no initial J and hence |3 = 0. Each run begins with 

a highly non- thermal equilibrium spectrum and the Fourier coefficients 

are then followed in time. The initial loading sets n(k) with 

]kl = 6 k . , 7 k . , 8 k . , 9k . , and 10 k . initially non-zero and 
1 • rmnr mur mm/ ^ mm-' mm 
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all other coefficients to zero. Tie magnitudes of the non-zero 
coefficients arc given, but the phases are chosen randomly. 

Run 2 uses ^2 wavenumbers (:‘i positive and 21 negative) with 
magnitudes which lie predominantly above k = 1. Run 5 consists of 
42 wove numbers with a significant lumber of them having magnitudes 
below k « 1. Run 7 has 540 wavenumbers with several in both wave- 
number ranges. 

The computing time necessary for a run varied from less than 
2 minutes for Run 2 to approximately 1 heur for Run 7 on the MFE 
CDCJ600 'computer. This drastic chcnge in computing time is not only 
because the system of equations is larger in Run 7 than in Run 2 but 
also because of the nature of the equations themselves. Equations 
(6)- — (5) are examples of stiff differential equations. 16,17 
Typically one has a set of stiff differential equations whenever the 
largest time scale in the problem is orders of magnitude larger than 
the smallest time scale in the problem. Good numerical techniques 
for stiff differential equations are one of the current frontiers' in 
numerical mathematics. Kumeiucal methods which perform excellently 
on non-stiff equations may fare poorly on stiff equations. Time steps 
in many standard schemes must he chosen small enough to fit every time 
scale in the problem, even if the faster scales are transient. The 
time step in our current problem was chosen to meet stability require- 
ments; the time step was halved whenever the errors in the conserved 
quantities indicated that numerical instabilities were beginning to 


arise. 
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Time-averaged electric field energy spectra were computed at 
various intervals during the course of a ran. Tlie final time-averaged 
spectra arc compared with the ensemble-averaged spectra in Eq. (51). 

Figures 1, 2, and' 3 show spectra for Funs 2, 3, and 7> 
respective] y. Each of the three figures consists of three panels: 
panel (a) shows the initial spectrum; panel (b) shows a time -averaged 
spectrum during the middle of 'the run, and pane'l (c) shows the final 
time-averaged spectrum. The solid curve in panel (c) of each figure 
is the theoretical prediction. In Figure 1, panel (h) shows the 
spectrum 'for Run 2 averaged between times t = 8,0 and t = 16.0 while 
the spectrum in panel (c) is averaged between times t = 24.0 and. 
t = 32.0. In Figure 2 the spectrum in panel (b) is averaged beWeen 
times t = 64.0 and t = 128.0 and in "panel (c) the spectrum is averaged 
between times t =,192.0 and t = 256.0 for Run 3- The times in Figure 
3 for Run 7 are t = 8.0 and t = l6.0 for panel (b) and t = 24.0 and 
t = 32.0 for panel (c). Every wavenumber is plotted in Figures 1 
and 2. Because the density of points in a logarithmic plot becomes 
so great at the high-k end of the spectrum for Run 7 only every other 
point is plotted for 3 < k < 10 and every fourth point is plotted 
for k > 10 in Figure 3 . 

Time histories of two typical Fourier coefficients are shorn in 

Figure 4 for Run 2. Figure 4(a) shows the real part of n(k = l), 

.or equivalently the imaginary part of E(k = 1), as a function of 

time. Figure 4(h) shows the real part of n(k =8), or equivalently 
* 

the imaginary part of 8E(k = 8), as a function of time. 
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Figure 5 shows time histories of two Fourier coefficients for 
Run 5. Figure 5(a) shows the real part of n(k = l/lo) vs. time and 
Figure 5(b) shows the real part of n(k = 1/2) plotted against time. 

In both Figures 4 and 5 the oscillations predicted by the 
linear theory seem to be present. A linear theory for Eqs. (6) — (?) 
gives, in our dimensionless units, the dispersion relation 

id 2 = 1 + k 2 (55 ) 

relating the frequency to the wavenumber k. For k = l/l6 this 
relation gives a period T = 2tt/«) =- 6.1. Counting peaks in Figure 5(a) 
and dividing into the total time interval between the first and last 
peak gives T ^ 6.2. A similar calculation for k - l/2 gives a 
theoretical period of T =- 5*6 and from Figure 5(b) T =* 5.8; for k = 1 
the theoretical period is T =* 4.4 and from Figure 4(a) T =* 4.1; for 
k = 8 the theoretical period is T =- O.785 while from Figure 4(b) 

T =- 0.778* It appears that the linear oscillations are present but 
that clearly something else is happening as well. 

Figure 6 shows the electric field energy as a function of time 
for the two modes of Run 2 shown in Figure 4. Figure 6(a) shows 

p O 

|E(k = l) [ versus time and Figure 6(b) shows jE(k = 8) j versus 'time. 
One should notice the orders of magnitude variation in the amplitudes 
of the modal energies. Because of this a time instantaneous wave- 
number spectrum has considerably more scatter in the points compared 
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to the tin. 0 -averc. 5 ed spectra of Figures 1 — 3. 
also shoe the more rapid flow of energy into rod out of llv.i h Lgher 
inodes- This can be expected from lis. ( 0 ) — ■('-•,) whore ve eon sec the 
the terms responsible for this energy transfer arc proportional to 
the wavenumber itself. Not only do the higher vavenunb-ar Fourier 
coefficients oscillate bn. a more rapid time scale but they also 
equilibrate faster as indicated by the spectre! plots in Figures 1‘ — 3 . 
This fast rush of energy to the large wavenumbers and hence smaller 
spatial scales would indicate the possibility for shock formation. 

Figures 7 , 8 , and 9 show configuration space plots of number 
density, velocity, and electric field, respectively, at four various 
times in ’the course of Run J. The shock formation is most clear in 
the velocity profiles of Figure 8 . Initially the velocity is zero. 

The first thing to occur is the linear transfer of energy from the 
electric and pressure fields to kinetic energy. At this point the 
nonlinear term in Eq. (7 ) begins to have effect. A simple analysis 
of the equation 


origwalpagb® 

os POOR QUAU™ 


Sv , &v . 
~ + V r- = 0 
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(34) 


shows that discontinuities form in a finite time. The time to this 
discontinuity given the data of Figure 8 at time t = 0.5 predicts the 
discontinuity to form at time t = 1 . 7 , later than the profiles indicate 
The fact that something drastic has happened is even more clear in the 



last density profile of Figure 7 in which the plot is virtually up 
and down between grid poJuLs. - Beyond time t = 1. 0 the configuration 
space profiles cannot mean much. This does not mean that the numeri- 
cal solutions of Eqs. (6) — (9) are incorrect or that the computed 
spectra cf Figures 1 — 3 are mathematically wrong. Because the time 
to discontinuity is much shorter than the time to equilibrium means 
only that the truncated Fourier representation in Eqs. (6) — (9) is 
no longer an adequate approximation to the continuous Eqs. (l) — (4) 
in the non-dissipative limit. We expect that with the addition of a 
finite dissipation at high wavenumbers that the discontinuities will 
not occur. If enough grid points are used (and hence more modes 
are kept in the truncation) we expect to be able to determine a shock 
width and the truncated Fourier equations will be adequate for very 
long times. The spectral laws may be different for the dissipative 
problem than they are here for the non-dissipative one just as the 
dissipative problem is significantly different from the non-dissipative 
one in Havier-Stokes fluids and MEED. What the above calculations do 
seem to indicate is that migrations of energy over large spatial 
scales is possible in a one-dimensional electron fluid. Any accurate 
assertions involving turbulent electron plasma oscillations must take 
these migrations into account. 



V. DISCUSSION 


In the previous sections we have produced a statistical inecfcan 
ical description} for turbulent motions in a one -dimensional electron 
fluid plasma. The prominent aspect of this paper i.s the application 
of the techniques of fluid turbulence- to a ’compressible situation. 
Much needs to be done., however., before applications to physical 
systems can he considered. The first obvious generalization is the 
addition 'of a finite dissipation to the problem. Dissipation will 
undoubtedly alter the results here as it has in -the analysis of other 
fluid equations. For example, Kolmogoroff arguments postulating an 
energy cascade would lead to jE(k) j ~ k~ • . Although migrations 
of energy across many length scales seems likely in this dissipative 
problem, it is unclear at this point whether the transfer of energy 
will be able to produce large scale electric fields through some 
sort of inverse cascade process, and there is no suggestion of them 
in the theory thus far. Shocks are also likely to play some role in 
this problem. 

Along other lines the inclusion of mobile ions could have 
unknown effects. Extensions to multiple space dimensions could also 
produce a variety of results. We already know that two-dimensional 
Navier-Stokes and MHD fluids behave differently from their three- 
dimensional counterparts. It is not inconceivable that a 
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multi-dimension electron fluid plasma may behave differently than 
the one-dimensional plasma studied here. 

We are also aware that Langmuir turbulence spectra which fall 

-2 

off as k at large khave arisen in simulations far more complex 

than this one (see, e.g. , Thomson et al.) and nay result from a 

« ^ 

variety, of arguments. We refrain from speculating as to whether the 
k behavior should characterize either (l) mobile ion situations, 
or (2) situations with finite dissipation. 
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APPET fulX A 

Frequently in studies of fluid equations via plane wave expan- 
sions the k = 0 mode is neglected. The k = 0 Fourier coefficient can 
he thought of as the mean value per unit length of the quantity under 
discussion. -For 'example 

E(k = 0) = jJ 1 E(x) dx 
n(k - 0) - ~ jJ 1 n(x) dx 

That the E(k = 0) mode cannot be automatically set to zero can be 
illustrated in the following example. Suppose we were to start our 
electron fluid initially with no number density fluctuations. Then 

E(x) = J X [1 - n(x)] dx 

implies that E(k = 0) = 0. Assume also that the initial velocity 
fluctuations are such to produce at a later time a surplus of elec- 
trons in one region and a deficit of electrons in another, main- 
taining overall charge neutrality. (See Figure 10. ) E(x) in the 
situation depicted in Figure 10 no longer has an average value of 



zero ami hence K(h = 0) ■/' 0. ll(h = 0) is just the jump in the 
scalar potential across one periodicity len£tii, and periodic electric 
fields do no l imply its vanishing. It is true that there are situa- 
tions in v.'lvich IS (It = 0) remains zero for all time, but these are not 
the most general situations one can envision. 





APPENDIX B 


In previous applications of truncated Fourier- series to Xavier - 
Stokes or KHD fluids the conserved quantities are all quadratic in 
the Fourier coefficients. In this appendix we will- demonstrate the 
conservation of S in Ea. (10) which involves products of three Fourier 
coefficients. There are four terms in £: the triples, the regaining 

kinetic energy, the mechanical energy associated with the pressure 
field, and the electric field energy. Vie trill differentiate each of 
the four terms, apply Eqs. (6) — (8) and sum the results. For this 
purpose it is helpful to rewrite Eas. (6) — (8) in the form 


= -ik £ 6 (p + r - k) n(p) v(r) - ikv(k) - ikv(o) n(k) 


pr^O 


(Bl) 


5v(k) 

at 


= -i % £ s(p + r - k) v(pj v(r) - ikv(O) v(k) 
Pi^O 


ikn(k) - ~ n(k) 


(B2) 


where ve have separated all terms with zero wavenumber coefficients 
and eliminated E(k) in favor of n(k) for k ^ 0. 6 (p + r - k) is the 
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Kronecker delta -which has the value 1 whenever its argument vanishes 
end the value 0 otherwise. Differentiating the triples we obtain 


1 c-. f Sv(k,) 


Sv(k 2 ) 3n(k.,) 

+ n^) v(k 3 ) — ^ + v(k 2 ) v(k 5 ) — — - 


Since the labels by which we refer to wavenumbers is immaterial, we 
notice that with the interchange of the labels kg and k^ in the second 
term that it is identical to the first term. Consequently, after 
applying Eqs- (Bl) and (B2) we obtain 


(triples) 


£ fiC^ + k p+kx) nf^Jvfkg) 


ik^ 

-o- L 6(p + r-k,) v(p) v(r) - ik,v(o). v(k, ) 

2 pr^O 5 5 3 


- ikj n(kj) - — n(kj) 



55 


+ \ . E h (k, + k + K ) v(k ) v(k, ) 

2 v v n .in i- d j d J 


k i k aV 0 


-ik E &(P + r -k,) n(p) v(r) 

. 1 pr^O ^ 


- ik 1 v(k 1 ) - ik 1 n{-k 1 ) v(O-) 


(B3) 


which consists of seven terms. In the first term we sum over k^. 

Again since the labels hy which we refer to wavenumbers is immaterial, 
after we have summed over we replace the labels p and r with the 
labels kj and k^. In the fifth term, which also involves four 
Fourier coefficients, after we have summed over we replace p and 
r with k^ and k^, respectively. The result for those two terms is 


i E 6 (k.. + k + k* +k. ) n(k ) v(k ) v(k=) v(k. ) i(k, + k. ) 

k 1 k 2 k 3 k 4 /0 


+ | E sC^ + kg + kj +k 4 ) nC^) v(k 2 ) v(k^) v(k 4 ) iC^ + k^) . 
- k l k 2 W° 

In the last term we may exchange the label k^ with the label k^. The 
resulting overall sum of these terms w T ith factors n(fc^) v(k^) v(k^ ) v(k 4 ) 
contains a factor (k^ + kg + k^ + k^) fiCk^ + kg + k^ + k^) which 



always vanishes. lienee the first and fifth terms cf ( Bi } cancel. 

In a similar way the second and seventh terms of (B3 ) cancel; their 


sum is 


1 

2 


o. 

k^kjf'O 


6 C- : i + kg + 1^) n^) v(k 2 ) v(l^) v(0) (h x 


OV ^ 

'■'* V 


Again exchange one for a k^ to obtain a factor (k^ + k ? + k^) 

6 (k^ + kg + k~) which always vanishes. The sixth terra, after cyclic 
permutations of wavenumbers k^ -» can be written as 


-I' £ 

“iW° 


5 (k 2 + k 2 + k 5 ) v(k 2 ) v(k 3 ) ik-jV^) 


| £ 6(kj_ + k 2 + kj) i(k x + k 2 + k^) vQ^) v(kg) v(k.) 

k l k 2 k ^° 


which also vanishes. Hence in summary 


3t 


(triples) = -i £ f sC^ + kg +k^) nCk^) v(kg) n(ls^)(k^ +—•) 
k l k 2 k 3^° ^ 


(B4) 


Differentiating the remaining kinetic energy terms of £ gives 


by a similar process 



ro ] h 



T,. v(-k) 

k^O 


nv(lr ) 



tf E e(p + r -k) v(p) v(r) 

P^o 


- ikv(o) v(k) - ikn(k) - i n(k) 


= -i 


kjtO 


n{k) v(-k) [k + 


(B5) 


Differentiating the pressure terms gives 



|n(k)| S 



Sn(k) 

St 



-ik E 6(p + r-k)n(p) v(r) 
. P^O 


- ikv(k) - ikv(o) n(k)l 


j 

- i E 6 (k x + k 2 + k 5 ) n(k^) v(kg) n(k^ ) k^ 

k^kgkj/o 

+ i E n (k) v(-k) k . (B6.) 
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.And finally the electric field energy terms give 



= £ ~~ 3 t^ -ih £3 6 (p + r - k) n(p) v(r). 

k~C 'K L piy'O 

- ikv(k) - ilcv(O) n(k) 

= i. 53 + nO^) v(J^) n(k 5 )~ 

k i k 2 k 3 ^° 5 



n(k) v{-k) | 


(B7) 


Summing E^s. (B4) — (B7) then gives Se/ftt = 0. ^j/dt = 0 can 


- be proven in a similar manner. 
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Table T 

List of Parameters for Three Typical 
Numerical Calculations 



Run 2 

Run 3 

Run 7 

N (number of grid points) 

64 

' 64 

512 

L (length of system in 2^) 

2tt 

32rr 

l6rr 

k . = Stt/l 

min ' 

1 

1/16 

1/8 

.k 

max 

21 

21/16 

170/8 

Size of time step 

( 64 )' 1 

(64 J' 1 

( 512 )" 1 
( 1024 ) _1 

Duration of run (in uT^*) 

pe 

32. 

256. 

32. 

Initial non- zero velocity 
coefficients 

Hone 

Hone 

Hone 

Initial non-zero density 
(electric field) 
coefficients 

k = 6, 7 , 
8, 9, 10 

k = 6/l6, 
7/16, 8/1 6, 
9/l6, 10/16 

k= 6/8, 7/8, 
8/8, 9/8, 
10/8 

Initial 6 
Pinal e 

Percent change in & 

0.134 
0.135 * 
6.67 XlO^ 

0.644 

0.644 

2.1x10 

0.260 
0.260 - 
8.2 xio" 6 

Initial J 
Pinal J 

0 _5 

-3.3 X10 5 

0 0 

-1.2X10 

0 g 

8.7 X 10" 

Initial Q 
Final Q 

0 -9 
1 . 5 X 10 * 

• ’ 0 -15 

4.3 x 10 

° -15 

-5.0 x 10 ? 

a 

0.00336 

0.0173 

0.0009 

Triples at end of run 

-o.oo 4 

-0.055 

-o.oi 4 
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FIGURE CAPTIONS 


Fig. 1 


Fig. 2 


Fig. 3 


Fig. 4 


Electric field energy spectra (E(k)j vs. k for Run 0; 

(a) initial conditions; (b) time averaged over tines t = 8.0 
to 16 . 0 ; (c) time averaged over times t = 24.0 to 32.0. 
Theoretical curve is the solid line in Fig. 1(c). 

Electric field energy spectra |e(1l) | vs. k for Sun 3: 

(a) initial conditions; (b) time averaged over times 
-t = 64.0 to 128 . 0 ; (c) time averaged over times t = 192.0 
to 256.0. Theoretical curve is the solid line in Fig. 2(c). 

p 

Electric field energy spectra |E(k) | vs. k for Sun 7; 

(a) initial conditions; (b) time averaged over times 
t * 8.0 to 16 . 0 ; (c) time averaged over times t = 24.0 
to 32.0. Theoretical curve is the solid line in Fig. 3(c). 
For clarity only every other value of k is plotted for 
3 < k < 10 and every fourth value for k > 10. 

Time histories of two typical Fourier coefficients in Run 2: 
(a) the real part of n(k = 1, t); (b) the real part of 
n(k = 8, t). 

Time histories of two typical Fourier coefficients in Run 3: 
(a) the real part of n(k = l/l6, t); (b) the real part of 
n(k = l/2, t). 


Fig. 5 



Fig. 6 Time histories of electric field energies for t\ o modes 

in Run 2: (a) jE(k - 1, t)| 2 , (b) |5(k = 0, t) | 2 , • 

Fig. 7 Number density profiles at four different t3>r.es in the 

evolution of Run 7* 

Fig. 8 Velocity field, profiles at four different tines in the 
evolution of Run J , 

Fig. 9 Electric field profiles at four different times in the 
evolution of Run 7- 

Fig. 10 Illustration showing the necessity for the E(k = 0) Fourier 
coefficient. The top graphs show the initial n(x) and E(x) 
with E(k = 0) s 0. The bottom graphs show the final n(x) 
and E(x} with E(k = 0) f 0. 
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Figure 8 













